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ABSTRACT 

A  new  reactive  flow  model  for  heterogeneous  energetic  materials  has  been 
developed  based  on  the  physical  and  chemical  parameters  of  the  material  as 
much  as  possible,  rather  than  solely  relying  on  empirical  constants  to  deter¬ 
mine  the  reaction  rates  behind  the  shock  wave.  Firstly,  this  report  presents  an 
extended  viscoplastic  pore  collapse  (hot  spot)  model  based  on  previous  models 
presented  in  the  literature.  Results  from  this  hot  spot  model  are  then  used 
to  develop  a  reactive  flow  model  embedded  into  the  multi-material  hydrocode. 
MULTI,  using  an  Induction-Parameter-Model  to  describe  the  thermomechan¬ 
ical  processes  of  pore  collapse  in  a  computationally  efficient  manner.  One 
and  two-dimensional  hydrocode  results  are  presented  for  the  energetic  ma¬ 
terial  HMX  undergoing  bare  and  cased  projectile  impact.  The  results  show 
the  importance  of  microstructure  in  determining  the  shock  ignition  and  sub¬ 
sequent  growth  behaviour  in  energetic  materials.  Thus,  a  new  capability  is 
described  for  determining  the  effects  of  varying  porosity  (due  to  manufacture, 
aging,  damage,  etc)  on  shock  sensitivity  and  can  be  used  to  help  evaluate  the 
Insensitive  Munitions  (IM)  qualities  of  weapon  systems. 

APPROVED  FOR  PUBLIC  RELEASE 


fo3-D$-/?§-3 


DSTO-TR-1383 


Published  by 

DSTO  Systems  Sciences  Laboratory 
PO  Box  1500 

Edinburgh,  South  Australia,  Australia  5111 

Telephone:  (08)  8259  5555 
Facsimile:  (08)  8259  6567 

©  Commonwealth  of  Australia  2003 
AR  012-547 
February,  2003 


APPROVED  FOR  PUBLIC  RELEASE 


u 


DSTO-TR-1383 


A  Microstructure  Dependent  Reactive  Flow  Model  for 
Heterogeneous  Energetic  Materials 


EXECUTIVE  SUMMARY 

A  reactive  flow  model  for  porous  energetic  materials  has  been  developed  using  a  visco¬ 
plastic  pore  collapse  model  to  describe  hot  spot  generation,  ignition  and  subsequent  burn¬ 
ing  through  a  pressure  dependent  burn  model.  The  viscoplastic  pore  collapse  model  is 
an  extension  of  models  from  the  literature  and  includes  an  enhanced  combustion  model 
with  a  BKW  equation  of  state  for  the  gases  trapped  within  the  pore.  Results  for  HMX 
indicate  that,  after  an  initial  thermal  explosion  within  the  pore,  a  quasi-steady  state  burn¬ 
ing  process  occurs  approximately  at  the  shock  pressure.  These  results  were  used  to  build 
a  reactive  flow  model  into  a  Eulerian  multi-material  hydrocode.  Instead  of  imbedding 
the  entire  viscoplastic  pore  collapse  model  into  the  hydrocode,  a  computationally  effi¬ 
cient  Induction-Parameter-Model  was  developed.  A  shock  initiation  experiment  for  HMX 
was  simulated  using  a  one-dimensional  model  and  it  was  found  that  a  tri-modal  initial 
pore  distribution  was  required  to  obtain  reasonable  agreement.  These  results  are  also  in 
agreement  with  another  reactive  flow  model  from  the  literature  that  solves  the  viscoplas¬ 
tic  pore  collapse  equations  at  every  time  step.  A  two-dimensional  simulation  of  a  single 
fragment  impact  experiment  was  also  performed.  These  results  illustrate  the  importance 
of  the  interaction  between  the  microstruture,  material  properties,  wave  structure,  cas¬ 
ing  and  fragment  impact  geometry.  Hence,  Weapons  Systems  Division  now  has  a  new 
modelling  capability,  enabling  shock  sensitivity  calculations  for  energetic  materials  with 
varying  amounts  of  porosity.  The  computationally  efficient  approach  used  here  lends  itself 
well  to  multi-dimensional  simulations  and  can  hopefully  be  used  for  Insensitive  Munitions 
(IM)  evaluations  of  future  and  in-service  weapon  systems. 
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1  Introduction 

This  report  documents  a  reactive  flow  model  for  energetic  materials  subjected  to  shock 
loading.  In  recent  years  shock  initiation  models  have  been  developed  which  are  relying 
less  on  empirical  data  to  calibrate  the  model  and  more  on  the  physical  and  chemical  prop¬ 
erties  of  the  material.  By  developing  such  models  shock  sensitivity  can  be  determined 
without  a  lengthy  experimental  program  to  calibrate  the  model  and  the  effect  of  changes 
in  microstructure  can  be  determined  computationally.  For  example,  the  effect  of  dam¬ 
age  (increase  in  porosity)  or  change  in  particle  size  can  be  determined  more  quickly  and 
inexpensively  than  by  using  a  purely  experimental  approach. 

It  is  of  great  advantage  to  those  involved  in  the  insensitive  munitions  community  to  be 
able  to  predict  the  high  speed  impact  (shock)  response  of  munitions  with  varying  degrees 
of  porosity  or  damage.  Hydrocodes  are  the  computational  tool  most  commonly  used 
to  model  shock  initiation  in  energetic  materials  and  the  ignition  process  and  subsequent 
chemical  kinetics  are  described  using  a  reactive  flow  model  imbedded  within  the  hydrocode. 
At  present,  the  most  reliable  reactive  flow  models  available  are  empirical  and  the  most 
common  and  successful  is  the  Lee  and  Tarver  (1980)  Ignition  and  Growth  model.  In  order 
to  use  this  model,  an  extensive  series  of  experiments  are  required  in  order  to  determine 
empirical  constants  for  the  model.  In  the  case  of  trying  to  determine  the  effect  of  increased 
porosity  or  damage,  an  expanded  series  of  experiments  are  required  for  each  microstructure 
which  potentially  could  be  very  expensive  and  time  consuming. 


1.1  Shock  Initiation  of  Energetic  Materials 

When  a  shock  wave  passes  through  an  energetic  material,  it  can  respond  in  different 
ways  depending  on  the  intensity  of  the  wave.  If  the  wrave  is  relatively  weak,  the  energy 
transferred  to  the  material  is  too  small  to  initiate  reaction  and  no  response  is  recorded 
other  than  mechanical  damage.  If  the  wave  strength  is  increased  to  a  moderate  level, 
chemical  reactions  are  initiated  but  are  quenched  due  to  insufficient  initial  energy  or  from 
release  waves  from  the  sample/ projectile  boundaries,  which  quickly  drop  the  pressure  and 
temperature  of  the  reaction  zone.  Increasing  the  wave  strength  beyond  a  shock  initiation 
threshold  value  results  in  the  shock- to-detonation  transition  (SDT).  In  this  case  the  energy 
contained  within  the  wave  is  sufficient  to  support  sustained  chemical  reactions  in  the 
compressed  region.  In  time,  the  pressure  builds  and  drives  a  detonation  wave  through  the 
material. 

When  early  researchers  calculated  the  bulk  temperature  in  the  compressed  zone  during 
the  SDT  process,  it  was  found  that  the  temperature  was  insufficient  to  cause  the  reactions 
observed  in  the  energetic  material.  It  was  concluded  that  energy  was  concentrated  into 
localised  zones  of  increased  temperature  known  as  ’hot  spots’.  The  chemical  reaction 
rates  are  increased  at  these  locations  and  if  the  conditions  are  right  will  result  in  sustained 
reaction  of  the  material.  Bowden  and  Yoffe  (1952)  were  the  first  to  suggest  the  existence  of 
these  hot  spots  and  subsequently  there  has  been  a  large  amount  of  effort  spent  investigating 
the  SDT  process  culminating  in  a  certain  consensus  of  opinion  on  the  processes  occurring 
during  shock  initiation  of  heterogeneous  energetic  material.  This  modern  view  of  shock 
initiation  is  summarised  in  Fig.  1  and  can  be  classified  into  four  stages.  The  first  stage 
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is  during  the  initial  shock  passage  into  the  material  where  energy  is  localised  into  hot 
spots  (hot  spot  generation).  If  the  size,  number  and  intensity  of  hot  spots  are  sufficient 
and  the  shock  pressure  is  maintained  for  a  sufficiently  long  period,  chemical  reactions  are 
sustained  and  burning  occurs  at  the  hot  spot  sites.  This  represents  the  second  stage  of  the 
process  (ignition).  As  the  hot  spots  are  at  discrete  locations  within  the  material,  burning 
initially  occurs  on  internal  surfaces  with  the  reaction  front  moving  outward  representing 
stage  three  (internal  burning).  When  a  reasonable  proportion  of  the  material  has  been 
consumed  (estimated  at  24-25%  by  volume),  the  burning  surface  switches  from  internal 
burning  to  external  burning  (stage  four,  external  burning).  This  occurs  when  either  hot 
spots  ’burn-through’  and  burning  is  established  around  each  individual  grain  or  the  solid 
material  becomes  sufficiently  weak  and  fractures  under  the  increasing  pressure  into  smaller 
burning  pieces  suspended  in  the  product  gases.  In  reality,  a  combination  of  these  processes 
is  likely  to  occur. 
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Figure  1:  A  schematic  illustrating  the  processes  occurring  during  shock  initiation  in  het¬ 
erogeneous  energetic  material.  The  diagram  shows  an  idealised  projectile  driving  a  shock 
wave  through  an  energetic  material  sample. 


1.2  Hot  Spot  Mechanisms 

There  has  been  much  debate  and  speculation  over  the  physical  mechanisms  that  cause 
hot  spots  during  shock  wave  interaction  with  energetic  materials.  It  is  however,  generally 
accepted  that  hot  spots  are  caused  by  the  interaction  of  a  shock  wave  with  defects  in  the 
microstructure  of  the  energetic  material  and  the  formation  process  is  controlled  by  the 
mechanical,  thermodynamic  and  chemical  properties  of  the  material. 

Many  hot  spot  mechanisms  have  been  proposed  and  despite  the  large  amount  of  effort 
spent  researching  these  phenomena,  a  single  process  has  not  been  identified  as  the  main 
energy  localisation  mechanism  during  shock  initiation.  A  comprehensive  list  of  hot  spot 
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mechanisms  will  not  be  presented  here  however  the  reader  is  referred  to  Bonnet  and  Butler 
(1996)  and  Mellor  et  ah  (1995)  for  a  comprehensive  list  of  possible  mechanisms.  Most  hot 
spot  mechanisms  for  shock  initiation  studies  involve  some  kind  of  viscoplastic  work  and 
this  is  regarded  as  an  efficient  means  of  using  shock  energy  to  create  high  temperature 
localised  zones. 

The  introduction  of  voids  into  an  energetic  material  increases  its  shock  sensitivity. 
This  well  known  fact  is  used  by  explosive  manufacturers  to  increase  the  detonability  of 
their  products.  As  mentioned  previously,  damaged  energetic  material  also  increases  shock 
sensitivity  as  additional  voids  are  introduced  via  the  damage  process.  Therefore,  it  is 
suggested  that  shock  induced  void  collapse  is  a  likely  hot  spot  generation  mechanism  for 
porous,  heterogeneous  and  damaged  energetic  materials.  Again,  a  number  of  void  collapse 
models  have  been  postulated.  A  successful  technique  that  incorporates  many  physical 
processes  is  the  viscoplastic  pore  collapse  model.  This  technique  is  based  upon  the  work 
of  Carroll  and  Holt  (1972)  who  investigated  the  interaction  of  shock  waves  and  porous, 
ductile  metal.  Khasainov  et  al.  (1981)  extended  this  work  to  energetic  materials  as  did 
Partom  (1981),  Frey  (1984),  Kang  et  al.  (1992),  Bonnet  and  Butler  (1996)  and  Massoni 
et  al.  (1999). 

The  viscoplastic  void  collapse  model  assumes  each  void  is  a  spherical,  gas  filled  pore. 
When  the  shock  passes  over  the  pore,  it  begins  to  collapse.  Heat  generation  is  caused 
by  the  compression  of  the  cavity  gas  and  viscoplastic  shearing  in  the  highly  stressed  zone 
of  material  near  the  pore  surface.  Thermal  conduction  into  the  bulk  of  the  material 
moderates  the  temperature.  Material  phase  changes  and  chemical  reactions  are  also  taken 
into  account  to  predict  the  onset  of  burning  or  explosion. 


1.3  Empirical  Reactive  Flow  Models 

Hydrocodes  commonly  use  empirical  reactive  flow  models  to  determine  the  shock  sen¬ 
sitivity  of  energetic  materials.  Probably  the  most  widely  used  model  is  the  Lee  and  Tarver 
(1980)  Ignition  and  Growth  model.  This  has  been  very  successful  in  modelling  the  shock 
initiation  of  several  explosives  and  propellants.  It  is  usually  expressed  as  a  three  term 
model  (Tarver  et  al.,  1996), 

dF  ,  / 0  \x 

—  =  /(l  F)b[y  1  aj  +GX(1  F)cFdpy  +  G2(l  F)eF9Pz  (1) 

where  F  is  the  fraction  of  energetic  material  reacted,  t  is  the  time,  r  is  the  material 
density,  pa  is  the  initial  material  density,  P  is  the  pressure  and  the  remaining  variables  in 
Eq.  1  are  empirical  constants.  The  first  term  describes  hot  spot  ignition  by  igniting  some 
of  the  material  relatively  quickly  but  limiting  it  to  a  small  proportion  of  the  total  solid 
( Figmax )•  The  second  term  represents  the  growth  of  reaction  from  the  hot  spot  sites  into 
the  material  and  describes  the  inward  and  outward  grain  burning  phenomena.  This  term 
is  limited  to  a  proportion  of  the  total  solid  (Fcimax)-  The  third  term  is  used  to  describe 
the  rapid  transition  to  detonation  observed  in  some  energetic  materials. 

The  Ignition  and  Growth  model  works  well  for  single  shock  waves  but  must  be  used 
with  care  for  multiple  shock  initiation  (Tarver  et  al.,  1995)  and  cannot  predict  the  effects 
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of  particle  size  or  porosity  without  empirically  classifying  each  material  with  differing  grain 
size  or  porosity. 

The  forest  fire  reaction  rate  model  (Mader,  1998)  is  also  widely  used.  This  model  uses 
empirical  ’pop-plots’  or  the  run-to-detonation  distance  and  single  curve  build  up  principle 
to  describe  the  shock  initiation  process.  Johnson,  Tang  and  Forest  (1985)  incorporated 
some  of  the  principles  of  the  forest  fire  model  into  their  JTF  reactive  flow  model.  This 
empirical  model  recognised  that  the  temperature  and  distribution  of  hot  spots  was  impor¬ 
tant  and  required  an  estimation  of  initial  hot  spot  temperature  and  how  it  develops  as  a 
function  of  shock  pressure. 

The  CPeX  code  (Jones  and  Kennedy,  1991)  uses  the  reactive  flow  model  of  Kirby  and 
Leiper  (1985)  to  study  the  detonation  processes  in  non-ideal  explosives  such  as  PBXW-115. 
This  model  accounts  for  hot  spot  ignition  and  subsequent  burning  through  pre-determined 
stages  defined  by  Gaussian  functions.  While  successful,  empirical  data  is  required  to 
describe  the  various  stages  of  reaction  in  terms  of  time  constants  and  mass  fractions  of 
various  components. 

Other  empirical  models  are  also  available  (e.g.  Anderson  et  al.,  (1981),  Partom  and 
Wackerle  (1989))  however  they  require  specific  information  calibrated  from  experiments 
that  do  not  easily  allow  for  the  adjustment  of  porosity  or  grain  size. 


1.4  Mechanistic  Reactive  Flow  Models 

Mechanistic  reactive  flow  models  attempt  to  incorporate  the  physics  of  hot  spot  for¬ 
mation  with  a  burn  model  (kinetics)  to  describe  shock  initiation  and  detonation.  The 
advantage  of  these  models  is  that  a  more  complete  description  of  the  shock  initiation  pro¬ 
cess  can  be  gained  thus  enabling  a  deeper  understanding  of  the  physics  occurring  at  the 
meso-scale.  Hence,  the  user  can  investigate  the  effect  of  varying  microstructural  properties 
on  shock  sensitivity  without  resorting  to  empirical  fits. 

Unfortunately,  modelling  all  the  processes  occurring  within  the  energetic  material  is 
very  computationally  expensive.  Baer  et  al.,  (1998)  attempted  this  by  using  a  high- 
resolution  shock  physics  code  to  model  mesoscale  shock  interactions  in  a  random  packed 
array  of  HMX  crystals  representing  a  small  piece  of  explosive.  Hot  spot  generation  was 
shown  to  occur  at  the  microstructural  defects  in  the  material  due  to  shock  focussing  and 
plastic  deformation.  However,  only  a  simple  burn  model  was  incorporated  due  to  the 
high  computational  cost.  To  simulate  a  1.2  mm  x  0.9  mm  x  0.9  mm  sample  at  68% 
theoretical  maximum  density  (TMD)  required  the  entire  Sandia  National  Laboratories 
TFLOP  (1012  floating  point  operations  per  second)  computer  using  4500  processors  with 
104  MB/processor.  It  is  obvious  that  to  use  this  method  to  model  larger  amounts  of  ener¬ 
getic  material  at  higher  TMD  and  with  improved  chemical  kinetics  will  require  significant 
advances  in  computational  performance. 

Similarly,  Conley  et  al.  (1998)  perform  microscale  simulations  of  PBX-9501  using  a 
2D  Eulerian  hydrocode  and  a  mesh  derived  directly  from  a  micrograph  of  a  PBX  section. 
Important  results  were  obtained  on  hot  spot  temperature  and  size,  however  a  reactive  flow 
model  was  not  included  presumably  because  of  the  high  computational  cost. 
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A  less  computationally  demanding  approach  is  to  use  a  continuum  based  mechanistic 
reactive  flow  model  imbedded  into  an  existing  hydrocode.  In  this  type  of  model,  the 
mechanisms  of  hot  spot  formation  and  subsequent  burning  are  assumed  in  advance  and 
applied  throughout  the  energetic  material.  Recent  examples  of  this  type  of  mechanistic 
reactive  flow  model  can  be  found  in  Cook  et  al.  (1998)  or  Bennett  (1998).  Making 
appropriate  assumptions  about  how  the  hot  spots  are  formed,  ignited  and  burnt  can  lead 
to  a  considerable  reduction  in  the  amounts  of  computational  resources  required.  The 
disadvantage  is  that  a  single  hot  spot  formation  and  ignition  theory  cannot  be  agreed 
upon,  which  may  raise  doubts  about  the  validity  of  the  model. 

As  mentioned  previously,  the  viscoplastic  pore  collapse  model  of  hot  spot  generation 
has  been  used  with  success  by  a  number  of  researchers.  Plastic  deformation  has  also  been 
identified  as  an  efficient  method  of  localising  shock  wave  energy.  Massoni  et  al.,  (1999) 
have  developed  a  mechanistic  reactive  flow  model  based  upon  viscoplastic  pore  collapse. 
The  model  was  used  to  successfully  describe  one-dimensional  shock  initiation  in  pressed 
HMX.  Despite  the  assumptions  used,  this  model  is  still  rather  computationally  expensive, 
requiring  approximately  150+  equations  to  be  solved  simultaneously.  This  effectively  limits 
the  model  to  one-dimensional  problems  using  the  current  computer  systems  in  Weapons 
Systems  Division,  DSTO.  This  report  will  outline  a  new  reactive  flow  model  based  upon 
the  work  of  Massoni  et  al.,  (1999)  which  is  less  computationally  expensive  and  will  allow 
two-dimensional  simulations  of  shock  initiation.  This  provides  a  useful  simulation  tool  for 
those  assessing  weapons  systems  for  unplanned  stimuli. 


2  Viscoplastic  Pore  Collapse  Model 

2.1  Model  Description 

The  viscoplastic  pore  collapse  model  of  Kang  et  al.  (1992)  has  been  used  to  describe  the 
processes  occurring  during  hot  spot  formation.  This  model  has  been  modified  and  extended 
by  Bonnett  and  Butler  (1996)  and  later  by  Massoni  et  al.  (1999).  The  model  presented 
below  represents  another  extension  of  the  Kang  et  al.  (1992)  model.  Here,  a  more  detailed 
description  of  the  heat  release  due  to  chemical  reactions  is  used  in  conjunction  with  the 
BKW  equations  of  state.  The  results  of  this  model  are  used  to  develop  the  reactive  flow 
model  presented  in  the  next  section. 

Figure  2  illustrates  the  mechanisms  assumed  during  viscoplastic  pore  collapse.  Fol¬ 
lowing  previous  work,  the  porosity  of  the  energetic  material  in  question  is  matched  in  the 
micromechanical  model  by  using  the  equation, 


where  (ft  is  the  porosity  and  a,  b  are  the  inner  and  out  radii  respectively  as  described 
in  Fig.  2.  Hence  the  average  pore  size  in  the  bulk  of  the  energetic  material  is  used  to 
determine  a  while  the  bulk  material  porosity  (ft  is  used  with  Eq.  2  to  determine  b. 
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on  Kang  et  al.  (1992). 


The  shock  is  assumed  to  pass  over  the  pore  instantaneously,  setting  up  a  hydrostatic 
pressure  field  in  the  material  surrounding  the  pore.  This  assumption  is  thought  to  be  valid 
for  small  pore  sizes  up  to  a  few  tens  of  micrometres.  For  larger  pore  sizes,  asymmetric 
collapse  of  the  pore  occurs,  resulting  in  a  complex  flow  pattern,  sometimes  with  material 
jetting  across  the  pore.  Therefore  this  model  is  limited  to  pore  sizes  of  a  few  tens  of  /xm. 

Shock  pressures  experienced  by  energetic  materials  during  shock  initiation  are  usually 
greater  than  1  GPa.  In  this  pressure  range,  the  mechanical  strength  of  the  material  is 
quickly  overwhelmed  and  the  pore  collapses  under  the  action  of  rapid  viscoplastic  flow. 
As  the  pore  collapses,  intense  heating  occurs  near  the  pore  surface  under  the  combined 
action  of  plastic  and  viscous  flow.  Eventually,  the  temperature  increases  to  a  level  where 
phase  change  occurs  at  the  pore  surface.  There  is  some  debate  over  whether  the  phase 
change  involves  melting  at  the  surface  with  chemical  reactions  occurring  in  the  condensed 
phases  (Bonnett  and  Butler,  1996)  or  sublimation  directly  into  the  gas  phase  with  the  bulk 
of  the  energy  release  occurring  in  the  gas  phase  (Kang  et  al.,  1992,  Massoni  et  al.,  1999). 
Here  we  assume  that  sublimation  occurs  at  the  surface  with  all  combustion  occurring  in 
the  gas  phase. 

The  motion  of  the  pore  surface  is  described  by  the  following  equation  derived  by 
Massoni  et  al.  (1999), 


dt 


1 _ J 

1  d>3  j  l 


psa{  1  03) 

•2 


Pg  Ps  4/xsJ(l  0) 


-  —  )  2Fssi* 

Ps  Pg 


ign(0  In 

^m[4(^  +  £)(1  03)  £(1  03)] }  > 


(3) 


where  £  =  a  ( m/ps ),  t  is  the  time,  m  is  the  mass  flow  rate  at  the  pore  surface,  ps  is 
the  density  of  the  solid,  Pg  is  the  pressure  of  the  pore  gas,  Ps  is  the  shock  pressure,  ps  is 
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the  viscosity  of  the  solid  phase  and  Ys  is  the  yield  strength  of  the  solid  phase.  Equation  3 
is  very  useful  as  it  applies  the  viscoplastic  resistance  with  the  gas  pressure  against  the 
shock  pressure  to  provide  an  equation  of  motion  for  the  pore  surface  with  a  correction  for 
mass  flux  due  to  sublimation. 

The  mass  flux  condition  at  the  surface  is  defined  by, 


\ps  (d  us)}r=a  =  [pg  (d  ug))r=a  =  m,  (4) 

where  us  and  ug  are  the  velocities  of  the  solid  and  gas  phases  respectively, 
he  solid  phase  velocity  is  controlled  through  the  following  equation  of  motion, 


us(r)  = 


(5) 


Some  of  the  heat  generated  from  viscoplastic  work  and  gas  compression  is  conducted 
through  the  pore  material.  The  transient  temperature  profile  through  the  pore  material 
is  calculated  by  solving  the  unsteady  one-dimensional  spherical  heat  equation, 


dTs  _  i, 

P,s  C.s  —  ks 


d2Ts  2  dTs 
dr2  r  dr 


(6) 


where  cs  is  the  specific  heat  of  the  solid  phase,  t  is  time,  ks  is  the  conductivity  of  the 
solid  phase,  Ts  is  the  temperature  of  the  solid  phase,  r  is  the  radial  position  and 


*,  =  12^(^)2  +  2u(N),  (7) 

is  the  viscoplastic  heating  source  term.  Equation  6  is  solved  at  each  time  step  using  an 
implicit  Gauss-Seidel  numerical  procedure  with  successive  over-relaxation.  An  adiabatic 
boundary  condition  is  applied  to  Eq.  6  for  the  outside  of  the  pore  material.  For  the  inner 
surface,  the  boundary  condition  is  set  to, 


kg 


(8) 


where  kg  is  the  conductivity  of  the  gas  phase.  Tg  is  the  pore  gas  temperature  and  T, 
is  the  interface  or  inner  pore  surface  temperature.  Equation  8  was  proposed  by  Massoni 
et  ah  (1999)  who  performed  a  systematic  analysis  to  determine  an  approximation  for 
the  thermal  profile  across  the  pore  and  hence  an  estimate  of  the  thermal  boundary  layer 
thickness.  Once  the  inner  surface  starts  to  sublime,  the  inner  boundary  condition  for  Eq.  8 
is  set  to  zero. 

The  energy  of  the  gas  within  the  pore  is  increased  by  the  compression  imparted  by  the 
collapsing  wall  of  the  pore.  It  is  also  modified  by  the  mass  addition  during  sublimation 
and  by  chemical  reaction.  The  specific  total  energy  of  the  pore  gas  is  given  by, 
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Ns 

£  =  e  +  x: 

t=l 


(AHf) 


Tr,i 


Mi 


fit 


(9) 


where  E  is  the  specific  total  gas  energy,  e  is  the  specific  internal  energy,  yAlIjj^  . 
is  the  molar  heat  of  formation  of  species  i  at  the  reference  temperature  Tr  =  0  K,  Mi  is 
the  molecular  weight  of  species  i,  fi  is  the  mass  fraction  of  species  i  and  Ns  is  the  total 
number  of  species  in  the  gas  mixture.  The  heat  of  formation  term  provides  the  mechanism 
for  heat  release  and  absorption  during  pore  collapse. 

The  specific  internal  energy  can  be  written  as  a  function  of  temperature  and  the  mass 
fraction  of  each  of  the  individual  species  in  the  mixture, 


e  =  £/(«>  (T).  (10) 

i= 1 

Following  Sichel  et  al.  (2002),  the  specific  internal  energy  can  be  written  in  terms  of 
the  McBride  and  Gordon  (1967)  polynomials, 

e  =  ^2  Rifi  +  (ai,i  o~)  +  -j±Tk\  ,  (11) 

i=i  l  k= 2  K  ) 

where  a*,  are  the  polynomial  constants  (McBride  and  Gordon,  1967)  and  R{  is  the  gas 
constant  for  species  i.  Equation  11  has  been  modified  to  include  the  term  a  from  the 
BKW  equation  of  state  (described  later  in  this  report).  This  is  done  to  include  real  gas 
effects  that  occur  at  the  high  shock  pressures  present  during  pore  collapse. 

Specific  internal  energy  is  modified  by  compression  from  the  pore  wall  or  from  con¬ 
duction  into  the  energetic  material.  The  rate  of  change  of  specific  internal  energy  can  be 
written  as, 

e  =  —  (5 p  A  .  (12) 

apg  V  a  ) 

It  is  assumed  that  the  conduction  term  (the  first  term  in  the  brackets)  is  zero  once 
sublimation  begins. 

Chemical  reactions  are  assumed  to  occur  in  the  gas  phase  only  and  to  follow  a  single 
step,  irreversible  reaction  involving  twelve  species  (i.e.  Ns  —  12), 


Ns 

<y.\Z\  — >  y ~^a.iZi, 

i- 2 


(13) 


8 


where  Zi  is  the  chemical  symbol  for  species  i  and  ai  is  the  stoichiometric  coefficient  for 
species  i.  In  this  report,  calculations  are  performed  using  HMX  energetic  material.  Table  1 
lists  the  stoichiometric  coefficients  for  the  reaction  used  here.  These  were  determined 
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Table  1:  Stoichiometric  coefficients  for  HMX  reaction. 


i 

Zi 

Oti 

1  (Reactant) 

HMXfgas) 

1 

2 

°2 

2,905e-5 

3 

n2 

3.986 

4 

h2o 

3.472 

5 

co2 

0.564 

6 

CO 

3.398 

7 

H2 

0.412 

8 

NO 

1.789e-3 

9 

no2 

3.729e-9 

10 

ch4 

3.780e-2 

11 

ch3 

3.802e-4 

12 

nh3 

2.62e-2 

by  performing  constant  volume  explosion  simulations  using  the  CHEETAH  (Fried  and 
Howard.  2001)  software  package. 

The  gas  combustion  is  assumed  to  follow  a  single  step  Arrhenius  law. 

ui  =  f\A  exp  RTs ,  (14) 

where  Co  is  the  rate  of  combustion  of  the  reactant  (HMX),  A  is  the  pre-exponential 
factor,  Ea  is  the  activation  energy  and  R  is  the  universal  gas  constant. 

Sublimation  is  assumed  to  occur  once  the  material  reaches  its  melting  temperature, 
which  follows  a  linear  relationship  with  shock  pressure  (Massoni  et  al.,  1999), 

Tm  =  Tm0  +  0PS,  (15) 

where  Tm  is  the  melting  temperature,  Tm q  is  the  melting  temperature  at  atmospheric 
pressure  and  6  is  a  constant.  Determining  the  sublimation  rate  is  difficult  as  it  represents 
the  first  step  (pyrolysis)  of  the  combustion  process  and  hence  is  controlled  by  the  kinetics  of 
decomposition.  As  these  rates  are  difficult  to  determine,  a  pressure-dependent  sublimation 
law  is  assumed, 


m  =  psApavPg  ,  (16) 

where  Ap  is  the  pore  surface  area  and  av ,  n  are  constants.  The  assumption  of  sub¬ 
limation  occurring  at  the  melt  temperature  has  been  chosen  as  a  part  of  a  convenient 
description  of  the  combustion  process  first  suggested  by  Massoni  et  al.  (1999). 

The  specific  internal  energy  is  updated  by  determining  the  contributions  from  me¬ 
chanical  compression  and  conduction  (Eq.  12)  and  changes  in  heat  of  formation.  The 
temperature  is  calculated  at  each  time  step  by  solving  Eq.  11  iteratively  using  a  Newton- 
Raphson  technique. 
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A  Becker-Kistiakowsky- Wilson  (BKW)  equation  of  state  is  used  to  describe  the  state 
of  the  pore  gas  (Mader,  1998), 


<7  =  §^  =  1  +  Xexp/?X,  (17) 

Ulg 

where  Vg  is  the  molar  gas  volume  and, 

*  =  w+f’  (18) 

where  a,/3,6,K  are  parameters.  Here  we  use  the  parameters  given  by  Vaullerin  and 
Espagnaq  (1998)  for  nitramines.  The  value,  k,  is  the  average  covolume, 

Ns 

k  =  Y,fiki>  (19) 

i=l 

where,  in  this  case,  fa  represents  the  mole  fraction  of  the  gaseous  species  i  and  fcj  is 
a  constant  covolume  characteristic  of  species  i.  Equation  19  is  summed  over  the  gaseous 
species  only. 


2.2  Results 

The  viscoplastic  pore  collapse  model  has  been  used  to  investigate  the  physical  processes 
occurring  during  shock  ignition.  Calculations  were  performed  using  cyclotetra-methylene- 
tetranitramine  (HMX)  as  the  energetic  material  with  varying  amounts  of  porosity  and 
initial  pore  radii.  Table  2  lists  the  parameters  and  constants  used  in  the  pore  collapse 
model  for  HMX.  All  calculations  in  this  section  were  conducted  using  a  shock  pressure 
of  2.2  GPa.  For  all  simulations,  the  pore  was  initially  filled  with  oxygen  at  atmospheric 
pressure  and  at  a  temperature  of  298  K. 

The  parameters  listed  in  Table  2  were  obtained  from  various  sources  in  the  literature 
(Massoni  et  al.,  1999,  Bonnett  and  Butler,  1996,  Kang  et  al.,  1992).  While  some  param¬ 
eters  have  little  uncertainty  associated  with  their  value  (e.g.  density),  other  parameters 
such  as  viscosity  and  yield  strength  have  not  been  established  at  high  shock  pressures.  The 
parameters  used  here  therefore  represent  the  current  ’best  estimates’  for  these  parame¬ 
ters,  however  experimental  work  remains  to  be  done  to  properly  quantify  these  values. 
The  sublimation  law  parameters  ( av ,  n)  were  obtained  from  pressurized  burning  test  data 
in  the  literature  (Atwood  et  ah,  1999,  Kubota  and  Sakamoto,  1989).  This  data  is  limited 
to  350  MPa,  which  is  significantly  lower  than  the  shock  pressures  experienced  here.  The 
burn  law  is  therefore  an  extrapolation  of  current  models  and  should  be  updated  once  high 
pressure  burn  data  is  obtained  or  an  alternative  approach  is  available.  The  use  of  burn 
law  data  for  the  sublimation  law  is  an  assumption  based  on  the  modelling  presented  by 
Kang  et  ah  (1992)  where  an  asymptotic  deflagration  law  was  used  to  describe  pore  sur¬ 
face  mass  flow.  As  this  approached  worked  well,  is  has  been  replicated  here.  However,  a 
simplified  empirical  pressure  dependent  law  is  used  instead  of  the  asymptotic  solution  to 
reduce  computational  overhead  and  maintain  reasonable  accuracy. 
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Table  2:  Parameters  used  for  HMX  in  viscoplastic  pore  collapse  model. 


Parameter 

Description 

Value 

Po 

Initial  solid  phase  density 

1905  kgm  3 

Cs 

Solid  phase  specific  heat 

1031  Jkg  XK  1 

ks 

Solid  phase  conductivity 

0.5016  Wm  XK  1 

Us 

Solid  phase  viscosity 

65  kgm  xs  1 

Ys 

Solid  phase  plastic  yield  strength 

200  MPa 

Tm 

Melting  temperature  at  STP 

548  K 

Pm 

Melting  point  pressure  dependence 

1.8  x  10  7  kPa  1 

dy 

Burn  law  parameter 

2.5  x  10  7 

n 

Burn  law  parameter 

0.66 

Tg 

Gas  phase  ratio  of  specific  heats 

1.4 

kg 

Gas  phase  conductivity 

0.0833  Wm  !K  1 

E/R 

Activation  temperature 

22000  K 

A 

Pre-exponential  factor 

1.93  x  1016s  1 

a 

BKW  EOS  parameter 

0.5 

P 

BKW  EOS  parameter 

0.16 

9 

BKW  EOS  parameter 

400 

K 

BKW  EOS  parameter 

10.90978 

Figure  3  shows  pore  collapse  results  for  HMX  with  an  initial  porosity  of  0.5%  and  an 
initial  pore  radius  of  2  mm.  The  figure  shows  the  pore  internal  radius  collapsing  with  time 
under  the  action  of  the  shock  pressure.  The  gas  pressure  remains  negligible  until  just  prior 
to  the  point  of  minimum  pore  radius.  A  rapid  rise  in  gas  pressure  is  then  observed  wdiich 
corresponds  to  the  rapid  consumption  of  sublimated  HMX  from  the  pore  wall.  After  the 
rapid  rise,  the  pressure  reaches  a  steady  state  corresponding  to  a  steady  burning  period. 

Figure  4  show's  the  pore  interface  and  gas  temperature  histories  for  the  same  test  case. 
The  interface  and  gas  temperatures  are  identical  in  the  early  portion  of  the  compression 
process.  Once  the  melting  point  is  reached  at  the  interface,  sublimation  begins  and  the 
interface  temperature  remains  constant.  The  gas  temperature  continues  to  rise  gradually 
for  a  short  time  and  then  experiences  a  rapid  rise  as  the  HMX  is  rapidly  consumed.  The 
temperature  readies  an  approximate  steady  level  after  0.225  ms. 

Figure  5  illustrates  the  reaction  progress  during  pore  collapse.  After  sublimation  be¬ 
gins,  HMX  begins  to  build  in  the  gas  phase  inside  the  pore.  As  the  pore  collapses,  the 
thermodynamic  conditions  within  the  pore  change  until  a  thermal  explosion  occurs.  This 
is  indicated  in  Fig.  5  by  the  rapid  decline  in  HMX  mass  fraction  and  the  rapid  rise  in 
nitrogen  (product)  mass  fraction.  Once  this  condition  has  been  reached,  the  mass  frac¬ 
tion  of  HMX  in  the  pore  remains  very  small,  indicating  that  the  reaction  is  controlled  by 
the  pressure  dependent  sublimation  process  at  the  surface  of  the  pore,  rather  than  the 
temperature  dependent  chemical  kinetics  within  the  pore. 

Figures  6-10  compare  viscoplastic  pore  collapse  results  for  HMX  with  varying  porosity 
{<f>  —  0.5  3%)  but  with  a  constant  initial  pore  radius  of  2  yum.  As  can  be  seen  from  the 

figures,  changing  the  porosity  while  maintaining  the  initial  pore  size  results  in  only  minor 
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Figure  3:  Viscoplastic  pore  collapse  results  for  HMX  -  internal  pore  radius  and  gas  pressure 
history.  Ps  =  2.2  GPa,  <j>  =  0.5%  and  a0  —  2.0  pm. 


Figure  p.  Viscoplastic  pore  collapse  results  for  HMX  -  interface  and  gas  temperature  his¬ 
tory.  Ps  =  2.2 GPa,  cf)  =  0.5%  and  a0  =  2.0 pm. 


changes  in  the  collapse  dynamics  and  thermodynamic  conditions. 

Figures  11-15  again  compare  results  from  the  viscoplastic  pore  collapse  model.  For 
these  calculations,  the  porosity  is  maintained  at  1.5%  and  the  initial  pore  radius  is  reduced 
from  2  to  0.5  pm.  The  effect  of  reducing  the  pore  radius  has  a  more  significant  effect  on  the 
collapse  dynamics.  The  maximum  explosion  temperature  (the  peaks  in  Fig.  13)  decreases 
as  the  pore  radius  is  decreased  however  the  time  of  thermal  explosion  increases.  The 
steady-state  burn  temperature  is  achieved  more  quickly  with  a  smaller  initial  pore  radius. 
The  peak  HMX  mass  fraction  in  the  pore  gas  increases  and  occurs  at  a  later  time  as  the 
initial  pore  radius  is  decreased. 

By  examining  Eq.  3  and  following  the  arguments  presented  by  Bonnett  and  Butler 
(1996),  the  viscous  term  resisting  the  pore  collapse  motion  increases  in  magnitude  as  the 
initial  pore  size  decreases.  Therefore  as  the  initial  pore  size  decreases,  the  collapse  rate  is 
decreased  and  the  overall  viscoplastic  and  gas  heating  mechanisms  are  reduced.  This  in 
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Time  (microseconds) 

Figure  5:  Viscoplastic  pore  collapse  results  for  HMX  -  reaction  history.  Ps  =  2.2  GPa, 
4>  —  0.5%  and  a0  =  2.0  pm. 


Time  (microseconds) 

Figure  6:  Comparison  of  pore  collapse  gas  pressure  results  for  HMX.  P3  —  2.2  GPa, 
<j>  =  0.5  3%  and  aQ  =  2.0  pm. 

turn  reduces  the  maximum  temperatures  and  extends  the  point  where  thermal  explosion 
occurs.  Interestingly,  the  time  wdiere  the  pore  radius  reaches  a  minimum  (Fig.  12}  is  only 
affected  slightly  by  the  change  in  initial  pore  radius.  As  the  time  where  the  peak  pressure 
occurs  (Fig.  11)  is  only  slightly  affected,  the  time  where  the  pore  begins  to  expand  remains 
approximately  the  same. 


3  Reactive  Flow  Model 

3.1  Model  Description 

The  results  obtained  from  the  viscoplastic  pore  collapse  model  indicate  that  after 
an  initial  transient  combustion  event  (thermal  explosion),  the  pressure  and  temperature 
approach  a  steady  state  almost  immediately  after  the  minimum  pore  radius  is  reached. 
The  constant  conditions  and  low  mass  fraction  of  HMX  in  the  pore  gas  indicate  that 
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Figure  7:  Comparison  of  pore  collapse  internal  radius  results  for  HMX.  Ps  =  2.2  GPa, 
(j>  =  0.5  3%  and  aD  —  2.0  pm. 


Figure  8:  Comparison  of  pore  collapse  gas  temperature  results  for  HMX.  Ps  —  2.2  GPa, 
<f>  =  0.5  3%  and  a0  -  2.0  pm. 

the  rate  of  combustion  is  controlled  by  the  pressure- dependent  sublimation  from  the  pore 
surface.  The  results  also  indicate  that  the  steady  state  pore  gas  pressure  is  approximately 
the  applied  shock  pressure. 

These  results  can  be  used  to  develop  a  reactive  flow  model  for  incorporation  into  a 
Eulerian  hydrocode.  The  one-dimensional  Euler  equations  written  in  conservative  form 
are, 


dp  dpu 
dt  dx 
dpu  d  pu 2  +  P) 
dt  dx 

dpE  d  ( puE  +  Pu) 


dt 


dx 


0 

0 

0. 


(20) 


These  equations  are  solved  using  an  existing  multi-material  Eulerian  hydrocode  de- 
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Figure  9:  Comparison  of  pore  collapse  HMX  mass  fraction  results.  Ps  =  2.2  GPa,  <f  = 
0.5  3%  and  a0  —  2.0  pm. 


Figure  10:  Comparison  of  pore  collapse  nitrogen  mass  fraction  results.  Ps  —  2.2  GPa, 
<f  =  0.5  3%  and  a0  =  2.0  pm. 

veloped  within  WSD  (Jones  et  al,,  1998).  This  two-dimensional  code  solves  the  Euler 
equations  using  the  flux  corrected  transport  algorithm  (Oran  and  Boris,  1987).  Various 
equations  of  state  can  be  used,  however,  for  the  work  presented  here,  the  JWL  (Lee  et 
al,  1968)  equation  of  state  is  used  for  reactive  materials  (separate  forms  for  un-reacted 
and  reacted  components)  and  the  Mie-Gruneisen  (Cheret,  1993)  equation  of  state  for  non- 
reactive  materials.  Interface  tracking  between  materials  (such  as  the  projectile/explosive 
Interface)  is  performed  using  Youngs  method  (Youngs,  1982). 

Assuming  pore  burning  occurs  at  a  pressure  equal  to  the  applied  shock  pressure,  the 
evolution  of  the  porosity  after  Ignition  can  be  described  using  the  following  equation 
(shown  in  one-dimensional  form), 


dp4>  ,  dp<pu  d<f 
dt  dx  P  dt ' 


(21) 


This  equation  converts  the  porosity  as  the  shock  wave  processes  the  material.  The 
right  hand  side  represents  burning  of  the  material  through  porosity  change.  This  term 
can  be  calculated  three  ways.  The  first  represents  the  case  where  there  is  no  burning.  This 
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Figure  11:  Comparison  of  pore  collapse  gas  pressure  results  for  HMX.  Ps  =  2.2  GPa, 
<j>  —  1.5%  and  a0  =  0.5  2.0  pm. 


Figure  12:  Comparison  of  pore  collapse  internal  radius  results  for  HMX.  Ps  =  2.2  GPa, 
<f>  =  1.5%  and  a0  =  0.5  2.0  pm. 


can  occur  either  if  the  initial  shock  pressure  is  too  low  or  in  the  period  before  ignition 
during  viscoplastic  pore  collapse.  This  first  case  is  represented  by, 

^  =  0  /  <  1.0  (no  ignition).  (22) 

at 

An  Induction  Parameter  Model  is  used  to  determine  when  to  start  the  internal  burning 
within  the  pores,  where  /  is  the  induction  parameter.  Induction  Parameter  Models  have 
been  used  previously  to  represent  the  chemical  induction  time  of  kinetically  complex, 
supersonic  reacting  mixtures  such  as  gaseous  detonation  waves  (e.g.  Sichel  et  ah,  2002). 
In  the  present  work,  the  induction  parameter  is  used  to  represent  the  time  when  a  burning 
pore  reaches  its  minimum  radius.  Following  the  viscoplastic  pore  collapse  results,  evolution 
of  the  porosity  from  this  time  onwards  is  controlled  by  sublimation  from  the  pore  surface. 
Hence,  the  induction  parameter  contains  the  material  response  as  well  as  the  chemical 
reaction  information. 
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Figure  13:  Comparison  of  pore  collapse  gas  temperature  results  for  HMX.  Ps  —  2.2  GPa, 
d>  =  1.5%  and  a0  =  0.5  2.0  pm. 


Figure  14:  Comparison  of  pore  collapse  HMX  mass  fraction  results.  Ps  =  2.2  GPa , 
<f>  =  1.5%  and  aQ  =  0.5  2.0  pm. 


Initially,  the  induction  parameter  is  set  to  zero  throughout  the  energetic  material  and 
is  convected  using  the  equation, 


dpf  dpfu  =  p 
dt  dx  t  5 


(23) 


where  r  is  the  induction  time  calculated  from  the  viscoplastic  pore  collapse  model.  It 
wrould  be  very  computationally  demanding  if  the  Induction  time  were  to  be  calculated  for 
each  cell,  at  every  time  step  using  the  full  viscoplastic  pore  collapse  model.  Fortunately, 
the  induction  time  can  be  reduced  to  an  analytical  expression,  determined  from  prior 
computational  runs  using  the  viscoplastic  pore  collapse  model. 

Figure  16  shows  the  induction  time  for  HMX  calculated  using  the  viscoplastic  collapse 
model.  The  model  was  used  to  calculate  the  induction  time  for  pore  radii  varying  between 
0.5  and  2.0  pm  and  initial  porosity  varying  between  0. 5-3.0%.  Figure  16  represents  the 
distillation  of  126  computational  runs.  The  symbols  represent  the  mean  of  the  data  where 
the  error  bars  are  one  standard  deviation  in  the  computed  results.  For  most  cases,  there 
is  only  a  small  variation  in  the  computed  results  at  each  shock  pressure,  reflecting  the 
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Figure  15:  Comparison  of  pore  collapse  nitrogen  mass  fraction  results.  Ps  =  2.2  GPa, 
4>  =  1.5%  and  aQ  =  0.5  2.0  pm. 

viscoplastic  pore  collapse  results  presented  in  the  previous  section.  This  variation,  denoted 
by  the  error  bars,  becomes  larger  at  small  shock  pressures  and  almost  insignificant  at  high 
shock  pressures.  In  any  case,  the  variation  is  reasonably  small  and  for  the  sake  of  simplicity 
of  the  computational  model,  a  power-law  curve  fit  was  applied  to  the  mean  data  (shown 
in  fig.  16).  This  results  in  an  induction  time  model  for  porous  HMX, 

r  =  0.2173P  °5925.  (24) 

Equation  24  can  be  used  in  Eq.  23  to  determine  the  induction  parameter  /.  In  some 
situations,  the  applied  shock  pressure  fails  to  ignite  the  energetic  material.  Figure  17 
shows  the  ignition  go/no-go  curve  as  calculated  by  the  viscoplastic  pore  collapse  model.  If 
the  pressure  falls  below  the  curve  in  Fig.  17  and  the  induction  parameter  is  less  than  unity, 
the  energetic  material  is  not  expected  to  ignite.  The  curve  illustrates  that  the  ignitability 
(not  detonability)  of  the  energetic  material  is  inversely  proportional  to  the  initial  pore 
radius  and  only  weakly  dependent  on  the  porosity  at  each  initial  pore  radius. 

Again,  a  power-law  curve  fit  was  applied  to  the  computed  results, 

P  =  0.9710ao°'8842.  (25) 

During  the  solution  of  Eq.  23,  Eq.  24  and  Eq.  25  are  solved  to  determine  the  source 
term  (right  hand  side  of  Eq.  23).  If  the  pressure  is  less  than  the  pressure  calculated  by 
Eq.  25,  then  the  shock  pressure  is  too  low  to  cause  ignition  and  the  source  term  is  set  to 
zero  (or  the  induction  time  is  set  to  infinity) .  Otherwise,  Eq.  24  is  used  to  determine  the 
induction  time  and  used  in  the  source  term  of  Eq.  23. 

Referring  to  Eq.  21,  the  right  hand  source  term  needs  to  be  evaluated  once  ignition 
has  begun  (i.e.  when  /  >  1.0).  Burning  occurs  initially  on  the  inside  of  the  pore  surface. 
After  a  certain  amount  of  time,  the  burning  changes  from  an  internal  surface  burn,  to 
an  external  surface  burn.  This  happens  when  the  energetic  material  dislocates  and  hot 
gas  escapes  from  the  pore  to  become  the  medium  in  which  individual  grains  are  burning. 
Shock  initiation  pressure  records  (Massoni  et  al.,  1999)  show  that  only  a  change  in  the 
burning  mode  can  explain  the  observed  curvature  of  the  pressure  response.  Since  detailed 
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Figure  16:  Induction  time  (t)  versus  Shock  Pressure  for  HMX,  <p  =  0.5  3.0%  and 

a0  =  0.5  2.0  pm.  Error  bars  indicate  one  standard  deviation  in  computed  results.  Red 

line  is  a  power  law  curve  fit. 

observation  of  when  the  change  in  burning  occurs  is  impossible,  it  is  assumed  that  it  occurs 
at  the  ’fluidisation  limit’  (Lee  and  Tarver,  1980,  Massoni  et  al.,  1999)  which  is  defined  as 
the  maximum  relative  volume  of  a  compact  bed  of  spheres  and  corresponds  to  /  =  0.24. 

There  are  two  different  source  terms  for  Eq.  21  depending  on  the  burning  mode.  These 
are  (for  /  >  1.0), 


NporeSpavPn  (f>  <  0.24  (internal  pore  burning) 

NpartSmavPn  >  0.24  (external  pore  burning),  (26) 

where  Npore  is  the  number  of  pores  per  unit  volume,  Npart  is  the  number  of  particles 
per  unit  volume,  Sp  and  Sm  are  the  pore  and  particle  surface  areas  respectively  and  av 
and  n  are  the  sublimation  parameters  used  in  the  viscoplastic  pore  collapse  model.  In 
order  to  determine  the  surface  areas,  the  current  pore  or  particle  radius  must  be  known. 
These  values  are  calculated  using  the  following  expressions, 


df 

dt 

df 

dt 


; 

y47T  Npore  J 


Rpart  — 


3(1  f) 


4  7T  N part 

where  a  is  the  pore  radius  and  Rpart  is  the  particle  radius. 


(27) 
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Figure  17:  Go/no-go  ignitability  curve  for  HMX  as  calculated  by  the  viscoplastic  pore 
collapse  model. 


The  pore  and  particle  concentrations  must  be  conserved  during  processing  by  the  shock 
wave.  Solving  the  following  conservation  equations  ensures  this, 


9Np0re  duNpore  _ 

dt  dx 

dNpart  ^  duNpart  _  q  (28) 

dt  dx 

Initially  the  pore  and  particle  concentrations  are  uniformly  set  throughout  the  energetic 
material.  These  initial  values  are  set  by  using  the  expressions  shown  in  Eq.  27  and  assuming 
an  initial  pore  radius,  porosity  and  particle  size.  The  initial  porosity  used  in  Eq.  27  is 
the  porosity  achieved  after  shock  compression  when  the  pore  radius  is  at  a  minimum. 
In  general,  this  porosity  is  approximately  one-hundredth  the  porosity  of  the  pre-shocked 
material. 

An  equation  of  state  is  now  required  to  close  the  system  of  equations.  As  mentioned 
previously,  the  Jones-Wilkins-Lee  (JWL)  equation  of  state  is  used  for  energetic  materials 
for  the  present  work.  The  JWL  equation  of  state  takes  the  following  form, 

P  =  Ae  R'v  +  Be  R>v +  (29) 

where  V  is  the  relative  volume,  T  is  the  temperature  and  A ,  B,  Ri,R.2,u)  (the  Gruneisen 
coefficient)  and  Cv  (the  average  heat  capacity)  are  constants.  There  are  two  sets  of  con¬ 
stants  for  the  JWL  equation  for  each  energetic  material,  one  for  the  unreacted  state  and 
one  for  the  reacted  state  (or  products).  Solving  the  JWL  for  a  reacting  mixture  therefore 
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requires  knowledge  of  the  mixture  fraction  and  a  mixing  law.  Typically,  the  mixture  frac¬ 
tion  is  obtained  using  the  Lee  and  Tarver  (1980)  Ignition  and  Growth  model  Eq.  1.  For 
the  current  reactive  flow  model,  the  porosity  is  the  assumed  indicator  of  volume  mixture 
fraction  of  reacted  material  (gaseous  reaction  products).  In  fact,  the  definition  of  porosity 
is  the  volume  fraction  of  the  gaseous  component. 

This  mixing  law  assumes  that  the  relative  volumes  of  the  gas  and  solid  phases  in  the 
energetic  material  can  be  described  using, 

(pvg  +  (1  (p)v8  =  vtou  (30) 

where  vg,  vs  and  vtot  represent  the  relative  volumes  of  the  gas,  solid  and  mixed  phases 
respectively.  The  system  is  then  assumed  to  be  in  pressure  equilibrium  and  an  iterative 
procedure  is  followed  where  the  relative  volumes  of  the  gas  and  solid  phases  are  varied  in 
accordance  with  Eq.  30  until  the  JYV  L  equation  of  state  (Eq.  29)  for  each  phase  computes 
the  same  pressure  within  a  specified  tolerance. 


3.2  One-Dimensional  Results 

Equations  21-30  were  incorporated  into  the  existing  Eulerian  hydrocode  (Jones  et  ah, 
1998)  and  simulations  were  performed  for  shock  initiation  of  an  energetic  composition  con¬ 
taining  HMX  with  1.7%  porosity.  Numerical  results  are  compared  with  the  experimental 
results  presented  by  Massoni  et  al.  (1999)  who  performed  shock  initiation  experiments  on 
a  HMX  sample  whose  porosity  and  microstructure  was  determined  prior  to  experimenta¬ 
tion.  Massoni  et  al.  (1999)  also  used  their  reactive  flow  model  to  simulate  this  experiment 
and  results  -were  very  favourable.  Here,  the  new  reactive  flow  model  presented  above  is 
used  to  simulate  the  same  experiments. 

The  description  of  the  experimental  setup  in  Massoni  et  al.  (1999)  is  poor  as  the 
only  details  given  were  a  column  of  HMX  with  1.7%  porosity  was  subjected  to  an  impact 
imparting  a  25  kbar  initial  shock  pressure.  In  order  to  simulate  this,  a  one  dimensional 
approach  was  adopted,  using  a  5  mm  thick  lexan  projectile  striking  a  25  mm  long  column 
of  HMX.  Figure  18  shows  a  schematic  of  the  numerical  model,  showing  the  locations  of 
the  Lagrangian  pressure  gauges  at  4,  7,  10  and  13  mm  respectively.  Lagrangian  pressure 
gauges  were  simulated  by  convecting  the  initial  positions  of  the  gauges  In  accordance  with 
the  Instantaneous  velocity  field.  The  actual  pressure  levels  were  obtained  by  extrapolating 
between  the  cell  values  either  side  of  the  convected  gauge  location. 

A  non-reactive  Mie-Grneisen  equation  of  state  (Meyers,  1994)  is  used  to  describe  the 
lexan  projectile  material.  Table  3  lists  the  parameters  used  for  the  equations  of  state  in 
the  numerical  model.  The  HMX  JWL  equation  of  state  was  assumed  to  be  the  same  as 
the  one  derived  by  Tarver  et  al.  (1993)  for  LX-10.  LX-10  is  an  energetic  composition 
comprising  of  94.5%  HMX  and  5.5%  Viton  Binder.  It  Is  pressed  to  1.6%  porosity. 

The  results  of  three  computational  runs  are  displayed  In  Fig.  19.  These  results  were 
calculated  using  mono-modal  and  tri-modal  initial  pore  distributions  within  the  material. 
In  order  to  calculate  the  initial  pore  number  density,  the  following  expressions  were  used, 
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Figure  18:  Schematic  of  one- dimensional  numerical  shock  initiation  model  for  HMX. 


Nporei 


1  yporep 

where  the  subscript  i  represents  the  particular  initial  pore  radius  (o0i),  in  is  the  fraction 
of  total  pores  with  initial  radius  a0i ,  Np  is  the  number  of  pore  distributions  considered, 
and  the  subscript  TReacted  JWL  represents  the  total  number  of  pores. 

Table  4  shows  the  initial  pore  distributions  used  for  the  computational  runs  displayed 
in  Fig.  19.  For  all  runs  in  Fig.  19,  the  lexan  projectile  impact  velocity  was  1.1  kms  1  and 
20  computational  cells  per  mm  were  used  to  ensure  a  converged  result.  A  time  step  of 
10  9s  was  also  used.  For  the  external,  particle  burn  model,  an  initial  particle  radius  of  35 
pm  was  used  for  all  runs  (Massoni  et  al.  1999). 


—  ViNporer 

3</>  ( 


4?r  VEfii ’ 


(31) 


Table  3:  Equation  of  State  Data  used  in  one- dimensional  numerical  simulations. 


Lexan  Projectile: 
Mie-Griineisen  EOS 

HMX: 

Unreacted  JWL 

HMX: 

Reacted  JWL 

p0  =  1.193  gem  a 

p0  =  1.865  gem  3 

c0  =  1.933  kms  1 

A  =  95220  GPa 

A  =  880.7  GPa 

s  =  2.04 

B  =  5.944  GPa 

B  =  18.36  GPa 

7o  =  0.61 

Ri  =  14.1 

Ri  =  4.62 

R2  =  1.41 

R2  =  1.32 

u  =  0.8867 

to  =  0.38 

Cv  =  2.7813  x  10  3  GpaK  1 

Cv  =  1.0  x  10  3  GpaK  1 

T0  =  298  K 

E0  =  10.4  GPacm3cm  3g 
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The  computational  results  indicate  that  the  shock  initiation  and  associated  pressure 
build-up  of  HMX  is  sensitive  to  the  initial  pore  size  distribution  within  the  material. 
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Table  4:  Initial  pore  distributions  for  computational  runs. 


Run  (a): 

Tri-modal 

Run  (b): 
Mono-modal 

Run  (c): 
Mono-modal 

aQl  —  0.5  pm  r)i  =  0.50 

a0l  =  0.35  pm  771  =  1,0 

aQ j  =  0.5  pm  r)i  =  1.0 

a02  =  0.35  pm  772  =  0.49 

aQ2  =  0  pm  772  =  0 

aD2  =  0  pm  772  =  0 

a03  =  0.2  pm  773  =  0.01 

a„3  =  0  pm  773  =  0 

a03  —  0  pm  773  =  0 

Npi  =  2.4276  x  1010  cm  3 

NP1  =  9.4658  x  1010  cm  3 

Nm  =  3.2468  x  1010  cm  3 

NP2  =  2.3791  x  1010  cm  3 

> 

13 

to 

II 

0 

£ 

to 

II 

0 

Np 3  =  4.8553  x  108  cm  3 

—  0 

CO 

II 

0 

Rpart0  —  35  pm 

Rpart0  —  35  pm 

Rparto  =  35  pm 

Fig.  19(a)  shows  reasonable  agreement  between  the  computational  model  and  the  experi¬ 
mental  results.  The  correct  shock  speeds  are  computed  along  with  the  initial  pressure  rises. 
The  peak  pressure  for  the  gauge  at  13  mm  is  under-predicted,  however  this  does  not  seem 
to  affect  the  run-to-detonation  distance  (shown  later).  When  the  simulation  is  repeated 
with  a  mono- modal  initial  pore  distribution  of  0.35  /xm  (Fig.  19(b)),  the  simulated  pressure 
build  up  occurs  much  more  quickly  with  an  associated  transition  to  detonation  observed. 
Similarly,  when  the  mono-modal  initial  pore  size  is  increased  to  0.5  /xm  (Fig.  19(c)),  a 
slower  pressure  build-up  occurs.  In  this  case,  excellent  comparisons  are  observed  for  the 
gauges  located  at  4  mm  and  7  mm.  however  reaction  rates  are  too  slow  to  maintain  the 
pressure  rises  for  the  gauges  at  10  mm  and  13  mm. 

These  results  reach  the  same  conclusions  as  the  study  by  Massoni  et  al.  (1999)  in  that 
to  model  the  shock  initiation  of  energetic  materials  accurately,  a  good  knowledge  of  the 
microstructure,  mechanical  and  chemical  properties  is  necessary.  In  the  study  by  Massoni 
et  al.  (1999),  a  tri-modal  initial  pore  radius  distribution  of  60%  0.5  pm,  39%  0.2  pm 
and  1%  0.15  pm  was  found  to  model  the  shock  initiation.  This  (small)  difference  in  pore 
distribution  can  probably  be  attributed  to  the  differences  in  the  viscoplastic  pore  collapse 
models  used  in  the  two  studies.  Here,  a  different  reaction  mechanism  and  combustion 
model  was  used  along  with  different  parameters  for  the  pressure  dependent  bum  law.  A 
BKW  equation  of  state  was  also  used  within  the  viscoplastic  pore  collapse  model  'whereas 
Massoni  et  al.  (1999)  use  the  H9  equation  of  state.  However,  despite  these  differences 
and  the  use  of  a  different  numerical  approach,  similar  results  were  obtained  albeit  with 
a  slightly  different  pore  distribution.  The  actual  pore  distribution  is  difficult  to  measure 
and  for  this  particular  case,  measurements  by  Massoni  et  al.  (1999)  suggest  a  majority  of 
pores  at  0.5  pm. 

As  the  details  of  the  shock  initiation  experiment  were  not  extensive,  the  method  of 
shock  generation  used  in  the  simulations  is  not  identical  to  the  method  used  in  the  exper¬ 
iment.  Therefore,  the  arrival  times  of  the  release  waves  from  the  free-end  of  the  projectile 
will  differ  between  computation  and  experiment.  The  complete  wave  system  is  shown  in 
Fig.  20  which  is  a  space-time  (x-t)  plot  of  equally  spaced  contours  of  pressure. 

The  interface  between  the  projectile  and  the  HMX  is  initially  at  0.5  cm.  From  this 
point,  two  initial  shock  waves  are  formed  in  the  HMX  and  projectile  respectively,  with  a 
region  of  uniform  pressure  between  the  shock  waves.  As  the  shock  in  the  HMX  processes 
the  material,  pores  are  collapsed  and  burning  begins,  resulting  in  an  increase  in  pressure 
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and  acceleration  of  shock  speed  (marked  as  the  Transition  region  on  Fig.  20).  Eventually 
the  shock  wave  transitions  to  a  constant  speed  detonation  wave.  Release  (or  expansion) 
waves  are  generated  when  the  initial  shock  in  the  projectile  is  reflected  from  the  free- 
end.  The  release  waves  are  transmitted  into  the  HMX  and  affect  the  pressure  levels 
after  the  transitioning  shock-to-detonation  wave.  As  the  detonation  wave  travels  faster 
than  the  release  waves,  the  run-to-detonation  distance  should  not  be  affected.  This  is 
illustrated  by  Fig.  21,  which  shows  a  comparison  between  the  run-to-detonation  distance 
(or  pop-plot)  measurements  of  Vanpoperynghe  et  al.  (1985)  and  the  numerical  model. 
Following  Massoni  et  al.  (1999),  the  run-to-detonation  distance  was  taken  to  be  the 
distance  along  the  explosive  where  the  slope  of  the  shock  wave  changes  on  a  space-time 
(x-t)  plot.  Initial  shock  pressures  were  varied  in  the  simulations  by  changing  the  impact 
speed  of  the  lexan  projectile.  The  experimental  measurements  are  represented  by  a  linear 
relationship  between  log  h  and  log  P.  The  energetic  composition  used  by  Vanpoperynghe 
et  al.  (1985)  was  96%  HMX  and  it  is  assumed  here  to  have  a  porosity  of  approximately 
1.7%.  The  comparison  between  experiment  and  calculation  is  reasonably  good. 

In  order  to  investigate  the  effects  of  porosity  on  shock  sensitivity,  two  simulations 
were  performed  using  the  tri-modal  pore  pore  distribution  used  in  Fig.  19(a)  but  with 
porosities  of  0.5%  and  3%.  As  expected,  the  porosity  level  has  a  significant  effect  on  shock 
sensitivity,  with  a  reduction  in  sensitivity  observed  for  the  lower  porosity  level.  A  similar 
change  in  sensitivity  can  also  be  observed  when  the  porosity  is  kept  constant,  but  the 
tri-modal  distribution  is  changed.  Figure  23  illustrates  this  effect,  showing  a  decrease  in 
sensitivity  when  the  average  initial  pore  size  is  increased  and  vice-versa  when  the  average 
pore  size  is  decreased.  This  can  be  explained  in  terms  of  the  available  surface  area  for 
combustion.  Provided  that  the  shock  pressure  is  high  enough  for  ignition  to  occur  at 
the  pore  locations  (see  Fig.  17),  the  smaller  pore  radii  allow  a  higher  surface  area  for 
combustion  to  occur  and  hence  an  increase  in  sensitivity  is  observed.  There  is  a  limit  to 
how  much  the  sensitivity  can  be  increased  by  reducing  the  pore  radii  and  this  is  given  by 
Fig.  17  using  the  viscoplastic  pore  collapse  model. 


3.3  Two-Dimensional  Results 

In  order  to  demonstrate  the  capabilities  of  the  reactive  flow  model  and  Eulerian  code 
the  single  fragment  impact  experiment  (as  described  by  Collignon  et  al.,  1992)  was  simu¬ 
lated  using  the  HMX  model  described  previously  in  this  report. 

Two-dimensional,  axisymmetric  simulations  were  performed  using  the  procedure  out¬ 
lined  above.  The  full  set  of  equations  used  for  the  axisymmetric  (cylindrical)  case  is, 


dp  dpu  1  drpv 

dt  dx  r  dr 

dpu  d  pu2  +  P)  1  drpuv 
dt  ^  dx  ^  r  dr 
dpv  dpuv  1  dr  pv 2  +  P ) 
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where  r  Is  the  radius  from  the  axis  and  v  is  the  radial  velocity.  The  source  terms  are 
as  described  earlier  in  this  report. 

A  schematic  representing  the  single  fragment  impact  test  is  shown  in  Fig.  24,  Here 
a  14.32  mm  diameter  steel  projectile  (representing  a  fragment)  travelling  at  847  ms  1 
strikes  a  76.2  mm  diameter  sample  of  HMX  within  a  6.35  mm  thick  aluminium  casing  and 
coverplate.  As  the  model  is  axisymmetric,  only  half  of  the  system  needs  to  be  represented. 
A  reflective  boundary  condition  is  imposed  on  the  centreline  to  simulate  the  other  half  of 
the  material.  Transmissive  boundary  conditions  are  imposed  on  other  boundaries  of  the 
model.  An  ordered  mesh  is  used  with  4  cells  per  mm  in  the  x-direction  and  r-directions 
in  conjunction  with  a  time-step  of  1.5  x  10  9  s.  A  Mie-Grneisen  equation  of  state  was 
used  for  the  steel  projectile  and  aluminium  casing  (parameters  from  Meyers,  1994)  and 
the  JWL  equation  of  state  shown  in  Table  3  was  used  for  the  energetic  material. 

In  the  experiments  performed  by  Collignon  et  al.  (1992),  various  HMX  based  formu¬ 
lations  were  tested.  For  the  purposes  of  simulation,  it  was  decided  that  the  formulation 
LX-14(N)  was  closest  to  the  HMX  model  derived  in  this  report.  This  formulation  consists 
of  95.5%  HMX  and  4.5%  estane  binder.  However,  no  information  on  the  microstructure 
(i.e.  porosity  level  and  pore  size  distribution)  is  available  therefore  the  trimodal  Initial 
pore  distribution  used  in  the  previous  section  is  used  for  this  simulation.  The  impact 
speed  of  847  ms  1  was  determined  experimentally  as  the  lowest  fragment  Impact  velocity 
that  caused  the  explosive  to  detonate. 

Figure  25  shows  results  from  the  two  dimensional  simulation  of  the  single  fragment 
impact  test.  In  each  figure,  contours  of  pressure  (top  half)  and  contours  of  fraction  of 
reacted  explosive  (lower  half)  are  shown  in  order  to  illustrate  the  processes  occurring  within 
the  energetic  material.  Initially  (Fig.  25(a)),  a  relatively  mild  shock  wave  is  imparted  to 
the  explosive  after  passing  through  the  aluminium  coverplate.  Reactions  are  initiated  as 
shown  by  the  well-spaced  contours  of  reacted  explosive.  Later  in  time  (Fig.  25(b)),  the 
shock  wave  strength  increases  due  to  the  quickening  reactions  behind  it.  Eventually  the 
shock  wave  transitions  into  a  detonation  wave  (Fig.  25(c)),  which  is  evident  by  the  defined 
high  pressure  detonation  front  and  the  concentrated  reaction  zone.  The  detonation  wave 
eventually  reflects  from  the  outer  aluminium  case  (Fig.  25(d)),  deforming  it  in  its  wake. 

Thus,  the  implementation  of  the  reactive  flow  model  Into  the  two-dimensional  hy¬ 
drocode  represents  a  significant,  new  capability  for  Weapons  Systems  Division.  The 
Induction-Parameter-Model  allows  a  computationally  efficient  approach  of  using  a  mech¬ 
anistic  reactive  flow  model  for  multi-dimensional  simulations.  Therefore  the  combined 
effects  of  microstructure,  casing  and  fragment  impact  geometry  can  be  explored  and  give 
meaningful  results  for  those  involved  in  Insensitive-Munitions  hazard  assessments. 
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4  Summary,  Conclusions  and  Future  Work 

A  reactive  flow  model  has  been  developed  based  upon  a  viscoplastic  pore  collapse  model 
for  hot  spot  generation.  The  viscoplastic  pore  collapse  model  represents  an  extension  of 
the  model  proposed  by  Kang  et  al.  (1992)  and  subsequently  developed  by  Bonnet  and 
Butler  (1996)  and  Massoni  et  al.  (1999).  Results  from  the  viscoplastic  pore  collapse  model 
indicate  that  pore  collapse  is  more  sensitive  to  changes  in  initial  pore  radius  rather  than 
bulk  porosity  levels.  The  model  also  shows  that  during  collapse,  ignition  begins  with  a 
thermal  explosion  inside  the  pore  that  quickly  develops  into  a  quasi-steady  burning  process 
that  occurs  approximately  at  the  shock  pressure. 

Using  this  information,  a  reactive  flow  model  was  developed  using  a  similar  framework 
to  Massoni  et  al.  (1999)  but  using  a  more  computationally  efficient  Induction-Parameter- 
Model  which  is  used  to  simulate  the  initial  shock  collapse  and  ignition  time  of  the  pores  in 
the  energetic  material.  The  Induction-Parameter-Model  therefore  represents  the  mechan¬ 
ical  and  reactive  properties  of  the  material  during  pore  collapse.  The  reactive  flow  model 
was  imbedded  in  an  existing  Eulerian  hydrocode  and  used  for  one  and  two-dimensional 
shock  initiation  simulations  in  HMX.  One-dimensional  simulations  are  directly  compared 
with  experimental  results  and  show  that  the  initial  distribution  of  pores  is  important  to 
accurately  simulate  the  shock  initiation  process.  There  is  also  reasonable  agreement  with 
the  run-to-detonation  distance  results  of  Vanpoperynghe  (1985). 

An  example  of  a  two-dimensional  simulation  is  presented  for  the  single  fragment  im¬ 
pact  test.  The  results  show  a  prompt  shock  detonation  initiation  occurs  for  the  particu¬ 
lar  projectile  diameter  and  impact  velocity  and  illustrates  the  hydrocode’s  capability  for 
determining  the  interaction  of  energetic  material  microstructure,  wave  structure,  casing 
material  and  fragment  impact  geometry.  Two-dimensional  simulations  are  made  possible 
on  the  current  WSD  computer  system  by  using  the  Induction-Parameter-Model  as  it  sig¬ 
nificantly  reduces  the  computational  cost  compared  with  the  full  implementation  of  the 
viscoplastic  pore  collapse  model. 

Improvements  can  be  made  to  the  reactive  flow  model.  The  pressure  dependent  burn 
model  is  extrapolated  well  beyond  the  pressure  range  of  the  experiments  it  was  determined 
from.  It  is  difficult  to  imagine  a  constant  pressure  burn  experiment  being  devised  in  the 
near  future  that  can  operate  at  detonation  pressures  and  give  the  parameters  required 
for  the  burn  model.  Instead,  a  better  approach  would  be  to  include  the  chemical  kinetics 
into  the  viscoplastic  pore  collapse  model  and  the  reactive  flow  model.  Tarver  et  al.  (1996) 
provide  a  three-step  chemical  decomposition  model  for  HMX  that  involves  reactions  in  the 
solid  and  gaseous  phases  (including  a  phase  change  step).  This  model  can  be  implemented 
to  improve  the  accuracy  of  the  models.  Further  improvements  include  a  more  rigorous 
treatment  of  the  phase  change  during  melting  and  obtaining  more  accurate  mechanical 
property  data  at  high  shock  pressures.  Also,  accurate  multiple  shock  simulations  are  not 
possible  using  the  present  formulation.  A  conductive  heat  loss  model  for  the  collapsing 
pores  needs  to  be  implemented  within  the  hydrocode  before  this  is  possible. 

A  new  gas  gun  facility  (Doolan,  2001)  is  nearing  completion  at  Weapons  Systems 
Division  that  will  provide  new  experimental  data  to  compare  with  the  reactive  flow  model. 
It  is  also  hoped  that  the  facility  will  be  able  to  provide  better  estimates  of  mechanical  data 
such  as  yield  strength  and  viscosity.  Incorporating  damage  and  aging  models  would  also 
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be  useful.  Such  models  could  predict  the  changes  in  porosity  and  physical  properties 
experienced  during  mechanical  insults  and  the  natural  aging  process  in  order  to  determine 
the  changes  in  shock  sensitivity  during  the  service-life  of  the  weapon  or  if  damaged  during 
combat  or  handling. 
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Figure  19:  Numerical  and  experimental  pressure  records  for  HMX  with  1.7%  porosity. 
Dots  represent  experimental  results  while  coloured  lines  indicate  numerical  results. 
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Figure  22:  Numerical  and  experimental  pressure  records  for  HMX  with  varying  poros¬ 
ity  and  a  constant  tri-modal  pore  distribution.  Dots  represent  experimental  results  while 
coloured  lines  indicate  numerical  results. 
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Figure  24:  Schematic  illustrating  single  fragment  impact  test. 
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